% queue_dynamics_delayed_heuristic(5,30)
% queue_dynamics_delayed_heuristic(5,60,18,'figure/exp0.mat')
function queue_dynamics_delayed_heuristic(provisioning_time,hold_duration,nsteps,file)
% randn('seed',12345)
% lambda = randn(1,400);
% %lambda(100:200)=4+lambda(100:200);
% %lambda=lambda+10;
% % randn('seed',16576576576);  k=cumsum(randn(200,1)*5 ); k= k-min(k); plot(k)
% lambda = lambda + sin((1:400)/20)*10+30;
 lambda= workload().get_workload('data/day42_per_min.txt', 1, 23*60)';
 lambda = [0 lambda];
 lambda(nsteps+1:end)=0;

mu=18;

 
a = 1;
b = [1/4 1/4 1/4 1/4];
lambda_smooth = filter(b,a,lambda);


r=1; target_q=10;
x0=20; 


% 
inst0 = [0;0;0];


%provisioning_time=5; %5
%hold_duration=1;   %60
% A=[0 0 0
%       1 0 0
%       0 1 0];
stations=provisioning_time+hold_duration;
A=[zeros(1,stations-1)
      eye(stations-1)];
A=[A zeros(stations,1)];   

% Ain = [1
%             0
%             0 ]; 
Ain=zeros(stations,1); 
Ain(1,1)=1; 

%  T=180; 
% % %m=mean(lambda(:,2:T+1));
%  Vartheta0=zeros(stations,1); 
%  Vartheta0(provisioning_time+1,1) = 5; 
% cvx_begin
%     variables X(1,T+1) Vartheta(stations,T+1) U(1,T)  cost_r
%     Vartheta(:,2:T+1) == A*Vartheta(:,1:T)+Ain*U(:,1:T);
%     X(:,2:T+1) == X(:,1:T)+lambda(:,2:T+1)-[zeros(1,provisioning_time) ones(1,hold_duration)]*Vartheta(:,1:T)  ; 
%     X(:,1) == x0;      
%     Vartheta(:,1)== Vartheta0;  % ones(stations,1); %[10; 4; 4]; 
%     Vartheta>0; 
%      X(:,2:T+1)>0;
%     cost_r==r*sum(U(:,1:T));
%     minimize (sum_square(target_q*ones(1,T)-X(:,2:T+1))...
%                         +cost_r)
% cvx_end
% Xopt=X;
% Uopt=U; 
% optcost = sum(sum_square(Xopt(:,2:T+1)-target_q))...
%                         +r*sum(Uopt)

% subplot(3,1,1); stairs(Uopt); title('Uopt'); 
% subplot(3,1,2); plot(Xopt); title('Q opt');
% subplot(3,1,3); plot(lambda_smooth); title('lambda'); 
% figure; plot(([zeros(1,provisioning_time) ones(1,hold_duration)]*Vartheta)'); 
% subplot(4,1,1); plot(Xopt); title('Xopt') 
% subplot(4,1,2); plot(Uopt'); title('Uopt') 
% subplot(4,1,3); plot(Xall); title('Xall') 
% subplot(4,1,4); plot(Uall'); title('Uall') 
%
% subplot(3,1,3); plot(lambda(2:end))

kk=0; 

%clear
% MPC
%for T=20 %
T = 60; % horizon
% nsteps =  18; % number of steps 
kappa=[]; 
%Vartheta0=zeros(stations,1);
varthetaPerInst = zeros(stations,11);      
varthetaPerInst(provisioning_time+1,1) = 5;  
% varthetaPerInstSum=zeros(1,11);         
% VarthetaPerInst0=
%Vartheta0(provisioning_time+1,1) = 5;  
n=1; m=11; 
x = x0;  
% vartheta=Vartheta0;

Xall = zeros(n,nsteps); Uall = zeros(m,nsteps); VarthetaAll=zeros(1,nsteps); 
kappaAll=zeros(2,nsteps);  
alpha_est= 0.75;
v_cost=10; 
lambda_hat=0.1576;
D=[zeros(1,provisioning_time) ones(1,hold_duration)];
dur = []; 

% a_in = eye(11);
map=[]; 
for i=1:100
    map(:,i)=ceil(i./prices_amz().get_compute_units); 
end  

for i = 1:nsteps
    vartheta = varthetaPerInst * prices_amz().get_compute_units(); 
    fprintf('.');  
    tic
    cvx_begin quiet 
        variables X(1,T+1) U(1,T)  Vartheta(stations,T+1)  cost_r  lmbda(1,T)
        % here I substituted lambda(:,2:T+1) with [mean(lambda(100:200)) mean(lambda(100:200))]
        X(:,1) == x;           
       lmbda == lambda_hat*ones(1,T) %lambda(:,1+i:T+i); %lambda_hat*ones(1,T)
        X(:,2:T+1) == X(:,1:T)+ lmbda -D*Vartheta(:,1:T)*mu;        
         Vartheta(:,1)== vartheta;  % ones(stations,1); %[10; 4; 4]; 
         Vartheta(:,2:T+1) == A*Vartheta(:,1:T)+Ain*U(:,1:T);
         Vartheta>0; 
         %X(:,2:T+1)>0;
         cost_r==r*sum(U(:,1:T));
        minimize (sum_square(target_q*ones(1,T)-X(:,2:T+1))...
                         +cost_r...
                         +v_cost *  sum_square(X(:,T+1)-target_q))
    cvx_end     
    dur(i) = toc;
    Xall(:,i) = x; % u = U(:,1); 
    x = max(x+ lambda(:,i) -D*vartheta*mu,0);    
    u_before = U(:,1);       
      
    % given u and bunch of rules, give me heuristic u
% %         kk=ceil(2.5./prices_amz().get_compute_units()).*prices_amz().get_hourly_price()
% %         find(kk=min(kk))
    u_inst = zeros(1,11);  
    
  
  u_before=round(u_before);  
  if (u_before>0)
        a__=sum(varthetaPerInst)';
        ind=find(a__>0);
        if (~isempty(ind))
            a__(ind)=1;
        end
        m=size(a__,1);
        b__=tril(ones(m,m),-1);
        b__=diag(a__)*b__; 
        % cost of hetrogenity and type 
        cost=(sum(b__*diag(map(:,u_before))))' +  diag(map(:,u_before))*prices_amz().get_hourly_price(); 
        cost(1)=Inf; %never use micro instance 
        ind=find(cost==min(cost)); 
        inst_type = ind(1);
        u_inst(ind(1)) =  map(ind(1),u_before); 
  end    
      
     varthetaPerInst= A* varthetaPerInst+ Ain* u_inst;   
     u_after = u_inst; % * prices_amz().get_compute_units; 
    %     for j=stations1:11
    %        ind=find(u_inst>0);
    %         u_inst=zeros(1,11); 
    %         u_inst(j) = map(j,u_before);  
    %        a_in=zeros(11,1);  ; 
    %        varthetaPerInst= varthetaPerInst*ones( + [u_inst zeros(1,stations-1)]; 
    %     end
    
    Uall(:,i) = u_after';  
%    vartheta = A*vartheta+Ain*u_after; %(u_' *prices_amz().get_compute_units());
 %   VarthetaAll(:,i) = D*vartheta; 
    lambda_hat= alpha_est*lambda_hat+(1-alpha_est)*lambda(:,i+1);
end


mpccost = sum(sum_square(Xall-target_q))...
                        +r*sum(Uall);
% disp(sprintf('T=%d',T)); disp(dur); 
disp(sprintf('mpccost=%d',mpccost));               
disp(sprintf('cvx_slvitr=%d',cvx_slvitr));
%end
%  plot(Uall)
% figure
% plot(Xall);

 % prices_amz().get_hourly_price()./prices_amz().get_compute_units()   


 %subplot(4,1,2);
 figure; stairs(Uall'); title('Uall'); 
 %subplot(4,1,3); stairs(Uopt); title('Uopt'); 
 % subplot(4,1,4); 
 figure; plot(lambda(1:nsteps)); title('lambda'); 
 
 % subplot(4,1,1);
 %plot(Xopt); hold on; 
 figure; plot(Xall,'--r');  title('Q opt');
 %axis([0 nsteps 0 20])
 save(file,'lambda','Uall','Xall'); 
% Annotate the point (-pi/4, sin(-pi/4))
%text(-pi/4,sin(-pi/4),'\leftarrow sin(-\pi\div4)',  'HorizontalAlignment','left')
 
% figure; plot(([zeros(1,provisioning_time) ones(1,hold_duration)]*Vartheta)'); 
% subplot(4,1,1); plot(Xopt); title('Xopt') 
% subplot(4,1,2); plot(Uopt'); title('Uopt') 
% subplot(4,1,3); plot(Xall); title('Xall') 
% subplot(4,1,4); plot(Uall'); title('Uall')  
% subplot(3,1,3); plot(lambda(2:end)) 

kk=0







